Experimental boson sampling in arbitrary integrated photonic circuits 
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Photons naturally solve the BosonSampling 
problem sample the outputs of a multi- 

photon experiment in a linear-optical interferom- 
eter. This is strongly believed to be hard to do on 
a classical computer [lj, and motivates the devel- 
opment of technologies that enable precise control 
of multi-photon interference in large interferom- 
eters Here we report multi-photon experi- 
ments in a 5-mode integrated interferometer. We 
use novel three-dimensional manufacturing tech- 
niques to achieve simultaneous control of 25 in- 
dependent parameters that describe an arbitrary 
interferometer. We characterize the chip using 
one- and two-photon experiments, and confirm 
the quantum mechanical predictions for three- 
photon interference. Scaled up versions of this 
setup are the most promising way to demonstrate 
the computational capability of quantum systems, 
and may have applications in high-precision mea- 
surements and quantum communication [5|. 

Large-scale quantum computers hold the promise of 
solving otherwise intractable computational problems 
such as factoring [f|. Despite all the experimental effort, 
this prospect is still far from feasible in all the proposed 
physical implementations Q. It is thus very important 
to establish intermediate experimental milestones for the 
field. One such example is the recent study, by Aaron- 
son and Arkhipov [l[ , of the computational complexity of 
simulating linear optical interferometers. It is well known 
that a linear-optical quantum computer, composed only 
of passive optical elements (such as beam splitters and 
phase shifters), becomes universal for quantum compu- 
tation if adaptive measurements are possible What 
was shown in [l[ is that such a device, even without 
adaptive measurements, produces an output that is hard 
to simulate classically. This suggests a feasible exper- 
iment to demonstrate the computational capabilities of 
quantum systems, consisting essentially in observing the 
multi-photon interference of Fock states in a sufficiently 
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large multimodc linear optical interferometer. 

More precisely, in [l[ it was shown that a linear opti- 
cal quantum process, consisting of (i) input of photons 
in a Fock state, (ii) unitary evolution implemented only 
via beam splitters and phase shifters and (iii) simulta- 
neous photon-counting measurement of all modes, can- 
not be efficiently simulated classically up to some rea- 
sonable complexity assumptions. This became known as 
the BosonSampling problem Efficient simulation, in 
this context, is understood as an efficient algorithm run 
on a classical computer that outputs simulated outcomes 
with a probability distribution which is close to the one 
obtained experimentally. At the core of the result is the 
fact that systems of noninteracting bosons evolve accord- 
ing to permanents of matrices [lOj. Suppose we have m 
bosonic modes in the initial Fock state 
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where Sj and a\ denote, respectively, the occupation 
number and creation operator for mode i, and |Q) de- 
notes the vacuum state. A linear optical evolution can 
be described by a m x m unitary transformation U on 
the space of creation operators, which induces a unitary 
transformation Up on the (exponentially larger) Fock 
space. The probability amplitude associated with input 
\S) and output |T) = \b\t<z...t m ) is given by 
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where Us,t is the matrix obtained by repeating su times 
the i th column of U, and tj times its j th row and 
per (A) denotes the permanent of matrix A [12], which is 
defined similarly to the determinant, but without nega- 
tive signs for odd permutations of matrix elements. 

Despite the similarity in the definitions, the perma- 
nent and the determinant are surprisingly different with 
respect to their computational complexity. While the 
determinant of a n x n matrix can be calculated in poly- 
nomial time, the permanent was proven to be compu- 
tationally #P-hard [T^], a class of intractable problems 
which includes the more well-known NP-hard problems. 
Despite having their dynamics ruled by permanents, non- 
interacting bosons cannot be directly used to calculate 
permanents efficiently, as a typical experiment will have 
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Fig. 1: Layout of multimode interferometers, (a) Re- 
alization of an arbitrary 5x5 mode transformation via a 
network of beam splitters with different transmissivities ti. 
The blue and red boxes stand for different phase shifters, (b) 
Implementation of the same scheme adopting integrated pho- 
tonics. 



an exponentially large number of outcomes, each pre- 
dicted by a hard-to-estimate, exponentially small proba- 
bility associated to a permanent by Eq. @. 

In [l[ it was estimated that a system of ~ 20 photons 
in m ~ 400 modes would already take noticeably long 
to simulate classically. At present, the most promising 
technology for achieving this regime involves inputting 
Fock states in multi-mode integrated photonic chips [2|- 
Given the theoretical hardncss-of-simulation 
results, experiments of this type serve as a milestone for 
quantum computation, on the way to more ambitious 
goals such as efficient integer factoring. 

In this Letter we report the experimental implementa- 
tion of a small instance of the Aaronson-Arkhipov pro- 
posal, using up to three photons interfering in a randomly 
chosen, 5-mode integrated photonic chip. We have made 
two important choices which provably make the quan- 
tum experiment harder to simulate classically [l[: we 
avoid any structure by choosing a random interferometer, 
and the interferometer has more modes than the number 
of input photons. Implementing this arbitrary interfer- 
ometer also serves as a stringent test of our novel man- 
ufacturing techniques, which crucially rely on a three- 
dimensional interferometer design. This allowed us to 
verify that non-interacting bosons evolve according to the 
permanent of matrices of size up to 3 x 3. 

Any m-mode linear interferometer can be decomposed 
in basic linear optical elements (phase shifters and beam 



splitters) using the decomposition given in [l8|. The gen- 
eral layout of these decompositions is depicted in Fig. 
[T}a for the case m = 5. It consists of a network of 
beam splitters with different transmissivities t{ (where 
tf is the photon's transmission probability) , interspersed 
by phase shifters restricted, without loss of generality, to 
the [0, 7r] range, as discussed in the Supplementary In- 
formation. Unfortunately, building large interferometers 
out of these discrete elements tends to result in mechani- 
cal instabilities which have prevented the demonstration 
of even a symmetric, 3-mode interferometer that pre- 
serves quantum coherence. A more promising approach 
to obtain stable multi-mode interferometers involves the 
fabrication of this network of linear optical element s by 
integrated optical waveguides in a glass chip [l9|, |20( ] . 
Waveguides are fabricated using the femtosecond laser 
micromachining technique [2lL l22j , which exploits non- 
linear absorption of focused femtosecond pulses to in- 
duce permanent and localized refractive index increase 
in transparent materials. Arbitrary three-dimensional 
circuits can be directly written by translating the sam- 
ple along the desired path, keeping the velocity constant 
with respect to the laser beam. This masklcss and single- 
step technique allows fast and cost-effective prototyping 
of new devices, enabling the implementation of three- 
dimensional layouts that are impossible to achieve with 
conventional lithography Q . 

In the integrated optics approach the role of beam 
splitters is performed by directional couplers, devices 
which bring two waveguides close together to redistribute 
the light propagating in them by evanescent field coupling 
[23[ • The integrated optics analogue of the discrete 
component layout, depicted in Fig. [T}a, is shown in Fig. 
[ljb, where one can appreciate the one-to-one correspon- 
dence between elements in the two approaches. Our main 
challenge in implementing the integrated layout of Fig. 
Q}b is to independently control each of the 10 transmis- 
sivities ti and 15 physically relevant [0,7r] phase shifts 
oti,!3i of an arbitrary, 5-mode chip. This is because in a 
typical optical circuit (Fig. [2ja) changes in the coupler 
geometry to modulate the transmissivity will change the 
optical path (and the phase shifts), and vice versa. 

The phase shifters are implemented by deforming the 
S-bcnt waveguides at the input of each directional coupler 
in order to stretch the optical path. The profile of the 
S-bends is deformed by a suitable coordinate transforma- 
tion (see Supplementary Information) that stretches the 
curve in a smooth fashion, to avoid adding waveguide 
losses, and does not modify the overall footprint of the 
S-bend, to avoid affecting the transmissivity of the sur- 
rounding couplers. Figurc[2]-b shows both an undeformed 
and a deformed S-bend. To determine the amount of 
phase shift that can be introduced by deforming the S- 
bend, we fabricated several Mach-Zchnder interferome- 
ters with increasingly larger deformations parameterized 
by the parameter d (see Supplementary information). 
This allowed us to calibrate our process (Fig. [U-b) and 
to verify that a phase shift of up to ti can be introduced 
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Fig. 2: Independent control of the phase shift and transmissivity at each directional coupler, (a) Controlled 
deformation of the S-bent waveguide at the input of each directional coupler and coupling geometry allow independent control 
over the phase shift and transmissivity. (b) The deformation is given by a non-linear coordinate transformation, which is 
function of a deformation coefficient d. The graph shows the undeformed S-bend together with a deformed one; the experimental 
dependence of the induced phase shift on the deformation parameter d at A = 806 nm is provided, compared to the expected 
one. (c) Control over the transmissivity of the directional coupler is performed by modulating the coupling coefficient; this is 
achieved by changing the waveguide spacing in the coupling region by rotating one arm of the directional coupler out of the 
main circuit plane. A calibration of the transmissivity dependence on the rotation angle at A = 806 nm is provided, compared 
to the theoretical expectation. 



without causing additional losses to the device. 

Achieving an independent control of the transmissiv- 
ity of each directional coupler is even more difficult. Two 
parameters that change the transmissivity are the inter- 
action length and the waveguide spacing. Changing ei- 
ther parameter induces, as a side effect, a variation in 
the optical path that could bring about a significant, 
and unwanted, phase shift. We overcome this limita- 
tion by using our three-dimensional design capability to 
rotate one arm of the directional coupler out of the main 
circuit plane, as depicted in Fig. [2]-c. This rotation 
is an effective way of modifying the waveguide distance 
in the coupling region (which changes the transmissivity 
between paths) without affecting the path lengths (and 
phase shifts). We found that rotation by a few degrees 
already enables us to span the full range of transmissivity 

(Fig. rjk). 

To choose which chip to fabricate, we sampled a uni- 
formly random 5x5 unitary and found its decomposition 
into directional couplers and phase shifters (see the Sup- 
plementary Information for more details). This chip was 
manufactured and used in single, two-, and three-photon 
experiments using the sources and detection apparatus 
described in Fig. [31 As a first step we characterized 
our 5-mode chip by injecting single photons in each in- 
put port i and measuring the probability P e (i,K) of 
detecting it in output mode K. The distribution prob- 
ability obtained experimentally is shown in Fig. 4-a, 



together with the theoretical prediction P^{i,K) of the 
sampled unitary U . To quantify the agreement between 
theory and experiment we calculated the similarity be- 
tween the two distributions, defined for two probabil- 
ity distributions p, q as S — (X)i \fPi1i) 2 ■ We obtained 
^exp t — 0.946 ± 0.005, which provides a first confirma- 
tion of the proper functioning of the device. Each output 
probability corresponds to the absolute value of one ma- 
trix element of U* . 

To obtain a complete characterization of the imple- 
mented interferometer we probed the device with pairs of 
photons. This was done by simultaneously injecting two 
single photons on all ten combinations of two different 
input modes For each input combination we esti- 

mated the ten output probabilities of the photons coming 
out in two distinct modes (K, L). In all, this corresponds 
to doing 100 Hong-Ou-Mandel, two-photon interference 
experiments (2~4j . The experimental distribution proba- 
bility Pe Xp (i, j, K, L) is reported in Fig. 4-b together with 
the theoretical distribution P^(i,j,K,L) expected from 
the sampled U l . Each theoretically predicted probability 
is obtained by calculating the permanent of a 2 x 2 sub- 
matrix of [/*. We observe a good agreement between the 
experimentally obtained probabilities and those given by 
the permanent formula (2), as evidenced by the similar- 
ity Sl xp t = 0.901 ± 0.027, thus confirming good control 
over the chip's fabrication parameters. 

It is possible to reconstruct the unitary matrix cor- 
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Fig. 3: Experimental setup for the characterization of the chip, (a) Schematic representation of the interferometer, 
realized by laser writing technique on a glass substrate, (b) One-photon, two-photon and three-photon states, generated 
by parametric down-conversion, are injected in the interferometer. The apparatus is composed of the elements: Avalanche 
Photodiode (APD), Interferential Filter (IF), Polarization Controller (PC), Polarising Beam Splitter (PBS), Parametric Down 
Conversion (PDC), Half Wave Plate (HWP). Spatial delay lines are adopted to synchronize the three photons, (c) Single-, two- 
and three-fold coincidence detection at the output ports of the chip is performed to reconstruct the probability of obtaining a 
given output state realization. 



responding to the multimodc interferometer using only 
the data corresponding to these one- and two-photon 
experiments |25|]. We have applied an adapted ver- 
sion of the algorithm described in (25|, obtaining a re- 
constructed unitary U r (see Supplementary Information 
for details). The similarity between the predictions 
of our reconstructed U r and our experimental data is 
Si r = 0.990 ± 0.005 (single photon experiments) and 
&exp,r = 0.977 ± 0.027 (two-photon experiments), which 
indicates a good characterization of the unitary imple- 
mented experimentally. 

We have also probed the chip's behavior in the multi- 
photon regime, by inputting three single photons into 
modes 1, 3 and 5 of our interferometer, and measuring 
the probability ratios of all events in which we find pho- 
tons exiting the chip in three different modes. In Fig.4-c 
we compare three distributions: the ideal distribution Pf 
obtained from 17'; the distribution Pf arising from our 
reconstructed U r and the one taking into consider- 
ation the partial indistinguishability p of the sources we 
used (for more details refer to the Supplementary Infor- 
mation). Fig.4-d shows a good agreement between the 
distribution P^ and our experimental results Pf xp as fur- 



ther confirmed by the similarity between these two dis- 
tributions S^ xp rp = 0.983 ±0.045. As these probabilities 
are proportional to permanents of 3 x 3 submatrices of 
the corresponding unitary, this is an experimental confir- 
mation of the permanent formula (2) in the three-photon, 
five-mode regime. 

We have experimentally confirmed that the permanent 
formula that governs the quantum mechanical behavior 
of non-interacting photons holds for up to three photons 
interfering in a randomly chosen, 5-mode interferometer. 
Scaling up experiments of this type would provide strong 
evidence of hard-to-simulate quantum behavior, even in 
the presence of noise [26| . This would require developing 
integrated multi-photon sources [27j |. detectors [28| and 
improving the manufacturing process to minimize losses 
and increase accuracy in the specification of each optical 
element. The capability of implementing arbitrary uni- 
tary transformations may find other applications, such 
as the discrete Fourier transforms required in the origi- 
nal KLM scheme for linear optical quantum computation 
[H, and basis-changing unitaries used in quantum state 
and process tomography [29l ] . 
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Fig. 4: Experimental results, (a) One-photon probability distribution: theoretical distribution (P t x (i, K j), experimental 
distribution (P exp (i, K)) and reconstructed distribution (P r (i,K)) . (b) Two-photon probability distribution: theoretical 
distribution (Pi(i, j, K, L)), experimental distribution {P^ xp {i, j, K, L)) and reconstructed distribution (P^(i, j, K, L)) . (c) 
Expected three-photon probability distribution for an input state (10101): theoretical distribution (Pf), reconstructed distri- 
bution (Pr) and reconstructed distribution with partial distinguishability (Pr, p )- (d) Three-photon probability distribution for 
an input state (10101): experimental distribution (Pe Xp ), reconstructed distribution with partial distinguishability (P^p). 
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Methods 

Femtosecond laser waveguide writing. The opti- 
cal circuit for the 5-modc interferometer has been fabri- 
cated by direct waveguide writing with the second har- 
monic (A = 515 nm) of a femtosecond Yb:KYW cavity- 
dumped laser oscillator. Our technique consists in focus- 
ing femtosecond laser pulses by a microscope objective 
(0.6 NA, 50 x ) in the volume of a transparent borosilicate 
glass (EAGLE 2000 - Corning). Under suitable irradia- 
tion conditions (300 fs pulse duration, 1 MHz repetition 
rate, 120 nJ energy), this creates a localized refractive 



index increase that can be exploited to write buried opti- 
cal waveguides by translating the sample with respect to 
the laser beam at uniform tangential velocity of 20 mm/s 
(Acrotcch FibcrGLIDE 3D air-bearing stages). Average 
depth of the fabricated devices under the glass surface is 
170 /im. The footprint of the 5-mode integrated circuit 
is 42 mm x 0.7 mm. 

Experimental apparatus. Four photons are pro- 
duced by parametric down conversion by pumping a 2mm 
long BBO crystal by a 392. 5nm wavelength pump field 
[30| . The four photons are generated at 785 nm, spec- 
trally filtered by 3nm interferential filters and coupled by 
single mode fibers. One of them acts as the trigger for 
coincidence detection, while the other three are coupled 
inside the chip after passing through different delay lines. 
The output modes are detected by using multimode fibers 
and single-photon avalanche photodiodes. Coincidences 
between different detectors allow us to reconstruct the 
probability of obtaining a given output state. 

Reconstructing the chip's unitary. In order to 
reconstruct the unitary which best approximates the ex- 
perimental data, we used an adaptation of the method 
reported in [25| . The method in pBT ] obtains a unitary 
that approximates the experimental data, in the form of 
the full one-photon outcome statistics (sixteen indepen- 
dent parameters, in our 5-mode chip), and a sufficient 
subset of the two-photon statistics (sixteen additional in- 
dependent parameters). By choosing different subsets of 
two-photon statistics used by the method, we obtained 
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25 different unitaries, each of which best fits the subset of 
the data used to obtain it. We compared the predictions 
of these 25 unitaries with the full data, picking the one 
with best agreement. This served as the starting point 



for a numerical search to maximize agreement with the 
experimental data, resulting in the reconstructed unitary 
U r reported in the Supplementary Information. 
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I. RANDOMLY SAMPLED AND 
RECONSTRUCTED UNITARIES 



The matrix U l was sampled from the uniform, Haar 
distribution over 5x5 unitary matrices: 
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where the global phase was fixed so as to make real. 
We then decomposed U* as a product of matrices that act 
non-trivially on two modes only, each representing a set 
of one beam splitter and two phase shifters in the range 
[0,7r], as described in [l| and depicted in Fig. 1 of the 
main text. The use of [0, it] phase shifters was an adap- 
tation of the decomposition of [3], which originally used 
phase shifters in the [0, 27r] range. This was done to limit 
the phase shift (and waveguide deformation) introduced 
by each element, as this may lead to losses. In this new 
decomposition, only one of the two phase shifters at each 
beam splitter input branch needs to introduce a non-zero 
phase shift. Table U reports the parameters obtained in 
this decomposition. 

Note that all experimental outcomes are invariant un- 
der multiplication of U by a phase shifter at each input 
and output port. For this reason we have set a\, fix , P5, j3% 
and /3io equal to zero in Table |H 

We have also found a unitary that fits well the single- 
and two-photon data, in part following the recently pro- 
posed method of @]. In Q, the authors show that the 
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unitary matrix corresponding to a physical interferome- 
ter can be reconstructed, up to a round of single phases 
at the input and output, using only single- and two- 
photon output probability distributions. The construc- 
tion requires choosing the matrix elements in one line 
and one column of the unitary matrix as reference points, 
and calculating the remaining elements by using a cor- 
responding subset of experimental data consisting of 16 
single-photon probabilities and 25 two-photon Hong-Ou- 
Mandcl visibilities. If the physical process is perfectly 
unitary, this freedom in the choice of the reference ele- 
ments is irrelevant, as all 25 possible choices produce the 
same matrix. 

However, the actual experiment has noise and the data 
cannot be perfectly described by one unitary matrix. 
This means that each of the 25 choices of reference el- 
ements, each of which corresponding to using one par- 
ticular subset of the data, produces a slightly different 
unitary. For our reconstruction, we have obtained these 
25 unitaries and calculated how well each of them repro- 
duces the whole data set. The measure of this quality 
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TABLE I: Transmissivities, ti, and phases, on and related 
to the layouts in Fig. 1 of the main text, that result from the 
decomposition of the sampled unitary matrix. 



was the x 2 between the whole set of single- and two- 
photon output predictions for the unitary with respect 



to the corresponding experimental data. We then chose 
the matrix which minimized this measure as a starting 
point for a numerical optimization. This optimization, 
a combination of brute force and gradient search meth- 
ods, slightly modified the best unitary obtained by the 
previous method to fine-tune the agreement with the ex- 
perimental data. In our reconstruction algorithm, the 
non-unitary indistinguishability (q ~ 0.95) between the 
two photons has been taken into account. 



Finally, since the reconstruction method of only ob- 
tains the unitary up to a round of arbitrary phases at the 
input and output modes, we have multiplied each row and 
column of the optimized unitary so as to obtain the high- 
est gate fidelity with U t , so that the matrices are easier 
to compare to each other. The reconstructed unitary U r 
is found to be: 



U r = 
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-0.474-0.147* 
-0.206- 0.319* 
0.265 + 0.125* 
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Notice that, once again, we have fixed the global 
phase so as to make U[ 1 real. The gate fidelity with 
respect to the sampled unitary U r was found to be 
F = |Tr([/*J7 rt )|/5 = 0.95. 




Fig. 1: Definition of the relevant parameters in an in- 
terferometer node. Scheme of a 3D directional coupler 
constituting the building block of the 5-mode interferometer, 
where we define all the relevant geometric parameters. 



II. MODELLING ARBITRARY PHASE SHIFTS 
AND TRANSMISSIVITIES 

Controlled phase shifts between subsequent directional 
couplers are implemented by varying the optical path, i.e. 
by stretching the geometrical length of the connection 
waveguide. This is perfomed by applying to the S-bent 
waveguide path the geometrical transformation described 
in the following. An undeformed S-bend is described by 
a sinusoidal function of the kind: 

h /2tt \ 

y = -2 cos \T x ) (3) 

where h and L define the extension of the curve in the two 
coordinates (Fig. [TJ. Such curve must be deformed and 
stretched in a smooth fashion, to avoid adding waveguide 
losses, and without modifying the overall footprint of the 
S-bend, which would otherwise affect the position of all 
the other couplers of the network. Hence, the deforma- 
tion is operated by the following coordinate transforma- 
tion: 

x = x + asm I —x J (4) 

where d defines the entity of the deformation. The re- 
sulting length of the deformed S-bend can be calculated 
by numerical integration. As discussed and shown in Fig. 
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2-b of the Main Text, the obtained phase shift was ex- 
perimentally calibrated with respect to the deformation 
parameter d. A good agreement is observed between the 
expected phase shift, calculated from the nominal length- 
ening of the S-bend, and the measured one. The root 
mean square deviation of the experimentally obtained 
phase shift with respect to the expected one is 0.25 rad, 
equivalent to an error in the path-length of about 20 nm. 
Independent control on the coupler transmissivity is ob- 
tained, on the other hand, by a three-dimensional rota- 
tion of one of the coupler's arm. Indeed, this enables to 
modify the distance between the two waveguides with- 
out altering the shape (and length) of the S-bend, which 
would determine an unwanted phase shift. To this aim, 
an accurate modelling of the coupling as a function of 
the rotation angle is needed. The transmissivity Q of a 
directional coupler, can be expressed by 

T = sin 2 (kZ) (5) 

where k is the coupling coefficient and Z is the interac- 
tion length (length of the region in which the two waveg- 
uides are close to each other) (Fig. Q}. The coefficient 
k decreases exponentially with the interaction distance 
s (spacing of the two waveguides within the interaction 
length) [J]: 

k = Koe _7 o (6) 

in which no and so are suitable constants. These con- 
stants have been experimentally measured by fabricat- 
ing directional couplers in the plane with different s and 
by fitting the theoretical curve ([6]) to the experimentally 
measured transmissions. The values arc k = 42 mm -1 
and So = 2.4 /xm. On the other hand, the distance be- 
tween the two waveguides depends on the rotation angle 
a of the arm of the directional coupler according to: 

s 2 = h 2 + (h + s min ) 2 - 2h(h + s min ) cos a (7) 

where h is defined above and s m ; n is the minimum dis- 
tance between the waveguides when the two arms are 
both in the same plane (a = 0). Putting together the 
Equations §E§ , ([6|) , §7$ one reads: 

T = sin 2 (e log KoL - * vc,1 " c ' 2COSQ ) (8) 



with C\ — h 2 + (h + s m i n ) and C2 = 2h(h + s m - m ) for 
brevity. Thus, inverting the formula ([8]) one can retrieve 
the angle a needed for achieving a specific value of T: 




Note that in Fig. 2-c of the Main Text the theoreti- 
cal curve shown in the graph is exactly Eq. ©, thus 
displaying the excellent agreement of the experimental 
points with the model here discussed. 



III. TWO-PHOTON AND THREE-PHOTON 
MEASUREMENTS 

A. Hong-Ou-Mandel visibilities 

The characterization process is based on performing 
single-photon and two-photon measurements. The latter 
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Fig. 2: Two-photon Hong-Ou-Mandel visibilities, (a) 

Experimental visibilities Ve Xp (i,j;K,L) and (b) predictions 
V^ p {i,j\K,L) with the reconstructed unitary U r obtained 
with two-photon measurements. 

consist on injecting all possible combinations of two pho- 
tons in two different input modes (i, j) and by measuring 
all possible Hong-Ou-Mandel interferences Q when the 
two photons exit in different output ports (K,L). This 
corresponds to the measurement of a 10 x 10 visibility 



4 



matrix, where the visibility is defined as: 



V 2 (i,j;K,L) 



P 2 ^(i,y,K,V)-P 2 {i, 3 ;K,L) 
P^(i,j;K,L) 



(10) 



Here, P 2 (i,j;K,L) is the two-photon probability for in- 
puts and output (K,L), and P 2 ' cl (i, j; K, L) is the 
two-photon probability when the two photons are delayed 
out of the interference region. The complete set of exper- 
imentally measured visibilities V 2 xp (i,j; K, L) is shown in 
Fig. [5]-a and is compared with the prediction obtained 
with the reconstructed unitary matrix U r of the inter- 
ferometer [Fig. H]-b]. The agreement between the two 
matrices is quantified by the similarity S t 
defined as: 



exp, rp , whi ch is 



S 2 = 1 - 

*-^exp,rp 



E«#j S^l IQm; K, L) - V 2 p {i, 3 ; K, L) 



2n 



where n = 100 is the number of measured visibilities. 
This parameter achieved a value S 2 xprp = 0.943 ± 0.003 
in our experiment, showing the quality of the reconstruc- 
tion process. 



B. Modeling photon distinguishability 

The three-photon state is conditionally prepared by ex- 
ploiting the second order process of a type-II parametric 



down-conversion source. The broadband feature of the 
pump is responsible for a partial distinguishability be- 
tween the generated photon pairs due to the presence of 
spectral correlations in the biphoton wave-function. This 
partial distinguishability introduces a reduction of the in- 
terference visibility in the three-photon experiment, and 
can be modeled by considering an input state of the form: 



g = r\l, 0,1, 0,1)(1, 0,1,0,1| + 
+ (l-r)|la,0,l 6 ,0,l a )(l o ,0,l 6 ,0,l o |. 



(12) 



Here, r is the parameter which takes into account the 
indistinguishability between the three photons, and the 
indexes i = a, b label that the photon on input port 3 
belongs to a different photon pair with respect to the 
photons on input ports 1 and 5. The parameter r is pro- 
portional to the overlap between the spectral functions of 
the three photons, and can be written as r = p 2 , where 
p is the indistinguishability of two photons belonging to 
different photon pairs. The latter has been measured di- 
rectly by performing a Hong-Ou-Mandel two-photon in- 
terference in a symmetric 50/50 beam-splitter, leading to 
p = 0.63 ±0.03. 
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